## Import data from .RData
load("data/cont_snotel.RData")

# Import seasonal snow zone polygons
SSZ <- st_read('data/MODIS_Snow_Zones/SSZ_0cc.shp')
## Reading layer `SSZ_0cc' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\SSZ_0cc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357223 ymin: 1163534 xmax: -437223.4 ymax: 3151034
## Projected CRS: NAD_1983_Albers
PSZ <- st_read('data/MODIS_Snow_Zones/PSZ_0cc.shp')
## Reading layer `PSZ_0cc' from data source 
##   `C:\Users\wreis\OneDrive - Colostate\MS Research\R_Git\SNOTEL\data\MODIS_Snow_Zones\PSZ_0cc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357223 ymin: 1459534 xmax: -437223.4 ymax: 3149702
## Projected CRS: NAD_1983_Albers
# Group sites by site_id
cont_snotel_site <- cont_snow_data %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise()
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# transform files to CRS:4326
snotel_sites <- st_as_sf(cont_snotel_site, coords = c('longitude', 'latitude'), crs = 4326) #fix before next run
PSZ_4326 <- st_transform(PSZ, crs = 4326)

# clip SNOTEL Sites to PSZ
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
PSZ_snotel_sites = st_intersection(PSZ_4326, snotel_sites)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
mapview(PSZ_snotel_sites) + 
  mapview(PSZ_4326)
# Filter all SNOTEL to PSZ SNOTEL Sites Only
PSZ_cont_snotel <- cont_snow_data[cont_snow_data$site_id %in% PSZ_snotel_sites$site_id,]

# PSZ SNOTEL Site, Calculate spring days with new SWE and new SWE depth 
PSZ_cont_snotel_spring <- PSZ_cont_snotel %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))

# SNOTEL Analysis, calculate the number and amount of positive SWE between 
cont_snotel_spring <- cont_snow_data %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"),yday = yday(date), year = year(date)) %>%
  mutate(newsnow = ifelse(snow_water_equivalent>lag(snow_water_equivalent), snow_water_equivalent-lag(snow_water_equivalent), 0)) %>%
  filter(state != "SD", year >= 1979, ifelse(leap_year(year) == FALSE, yday >= 60 & yday <= 243, yday >= 61 & yday <= 244))

# PSZ determine average number of days with increased SWE per spring at each site
PSZ_cont_snotel_site <- PSZ_cont_snotel_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years, 
            newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Determine average number of days with increased SWE per spring at each site
cont_snotel_site <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), years = length(unique(year)), avg_spring_days = newsnow_days/years, 
            newsnow_depth = sum(newsnow, na.rm = TRUE), avg_spring_snow = newsnow_depth/years)
## `summarise()` has grouped output by 'site_id', 'site_name', 'state',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# View data
mapview(cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
        layer.name = "AVG Sping Days (All Sites)" ,crs = 4326, grid = FALSE) + 
  mapview(PSZ_cont_snotel_site, ycol = "latitude", xcol = "longitude", zcol = "avg_spring_days",
        layer.name = "AVG Sping Days (PSZ Sites)" , crs = 4326, grid = FALSE)
# Days of new snow at SNOTEL Sites in Western US
#Summarize the number of new SWE days and SWE depth at all sites per year
cont_snotel_site_yr <- cont_snotel_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a box plot for all sites
spring_snow_days <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_days)
#Summarize the number of new SWE days and SWE depth at PSZ sites per year
PSZ_cont_snotel_site_yr <- PSZ_cont_snotel_spring %>%
  group_by(site_id, site_name, state, year, latitude, longitude, elev) %>%
  summarise(newsnow_days = sum(newsnow>0, na.rm = TRUE), newsnow_depth = sum(newsnow, na.rm = TRUE))
## `summarise()` has grouped output by 'site_id', 'site_name', 'state', 'year',
## 'latitude', 'longitude'. You can override using the `.groups` argument.
# Create a box plot for all sites
PSZ_spring_snow_days <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_days, color = state)) +
  geom_boxplot()+
  ggtitle("Days of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_days)
# Spring Snow depth
spring_snow_depth <- ggplot(cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
  geom_boxplot()+
  ggtitle("Depth of increased SWE at all sites (1978-2022)")
ggplotly(spring_snow_depth)
# Spring Snow depth
PSZ_spring_snow_depth <- ggplot(PSZ_cont_snotel_site_yr, aes(x = state, y= newsnow_depth, color = state)) +
  geom_boxplot()+
  ggtitle("Depth of increased SWE at PSZ sites (1978-2022)")
ggplotly(PSZ_spring_snow_depth)
# spring_snow_days_yr <- ggplot(cont_snotel_site_yr, aes(x = year, y = newsnow_days, group = state, color = state)) + 
#   geom_line()
# 
# ggplotly(spring_snow_days_yr)